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ABSTRACT 

Exact and approximate analytical formulas are derived for the internal struc- 
ture and global parameters of the spherical non-rotating quasi-incompressible 
planet. The planet is modeled by a polytrope with a small polytropic index 
n < 1, and solutions of the relevant differential equations are obtained analyti- 
cally, to the second order of n. Some solved and unsolved problems of polytropic 
models are discussed as well. 

Subject headings: planets: internal structure — quasi-incompressible model — polytropes 

1. Introduction 

Many classical problems of the theoretical astrophysics: internal structure, effects of 
rotation, tidal effects, pulsations and stability of the stars/planets are solved mainly in the 
case of incompressible (or homogeneous) liquid. 

The point is that even the simplest (and basic) differential equation of internal structure 
- Lane-Emden equation (hereafter LEE) of index n is solved analytically only for three 
values of n = 0, 1,5, each of them having their own deficiencies - case n = corresponds 
to incompressible liquid, polytropic model with n = 1 has constant radius independent of 
its mass, and the star with n = 5 (and with n > 5) has infinite radius. Some 25 years ago 
SK ( Seidov & Kuzakhmedov, 1977 (SK77), 1978 (SK78)) had presented the new analytical 
solutions of the LEE for index n only slightly differing from 0, 1, and 5, and also the series 
form of solution of LEE of arbitrary index. Several authors used and extended the approach 
of SK for the non-rotating and slowly rotating models - Seidov (1978a,b, 1979a,b); Seidov, 
Sharma & Kuzakhmedov (1979); Mohan & Al-Bayaty (1980); Jabbar (1984); Caimmi (1987); 
Horedt (1987, 1990); Medvedev & Rybicki (2001). 

In this note I return to SK and present some analytical results for the quasi-incompressible 
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model of planet. The idea is to use the perturbation theory of differential equations: if we 
have the analytical solution of the ordinary differential equation (ODE) for some particular 
value of parameter n = no (say n = 0) then we may try to look for the analytical solution 
of the same ODE with n ~ n (in our case n ~ 0). 



2. Basic equation 

The basic equation is LEE of index n : 



i d / 2 d y\ _ 



r dx \ dx J 
with initial conditions y(0) = 1, y'(0) = 0. 

We look for solution y(x) in an interval from x = to the first zero x = X such that |/(X) = 
(hereafter to be concise I use shorthand s = VE). 

Three classical analytical solutions of eq. (f) are long known (see e.g. (Chandrasekhar 
1957)): 

1 

n = 0, y = 1- -x , X = s, fi = 2s, p c fp m = l; (2) 
o 

n = l, 2/=:^^ X = tt, ^ = 7r, p c /p m = ?r 2 /3; (3) 

71 = 5, y = (1 + ^x 2 y 1/2 , X -> oo, /i = v / 3, p c /p m -> oo. (4) 

In these equations, p = — X 2 y'(X), p c / p m = X 3 /3p; X, p are dimensionless radius and 
mass, and p c / p m is the central-to-mean density ratio. 



3. The perturbation method 

Consider eq. (1) as ODE depending on parameter n, then assuming n as a small 
parameter, n <^ 1, we expand the r.s. of eq. (1) to the second order of n: 



y = y + nyi + n y 2 , y n = 1 + n In (y ) + n + - In (y ) ) • (5) 

From eqs. (1, 5) we have three coupled ODEs for three functions y , y± : y 2 : 

1 d ( 2 dyo 
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]_d( 2 djn\ = y , = 

r dx \ ax J 

1 d ( 2 dy 2 \ f yi , 1 . 2 



Initial conditions in eqs. (6, 7, 8) are defined by the form of series expansion of the solution 
of LEE of arbitrary n at x = 0: 

1 2 n , n(8ra-5) fi n (122 n 2 - 183 n + 70) 8 

» = 1 -6^ + 5i^ + -w J ^ + - iP9! (9) 

Expanding eq. (9) to the second order of n, we have the series expansions for functions yi, y 2 
at x = 0: 

Note that > and y 2 < 0. 

Before solving eqs. (7, 8) I'd like to dwell some on the validity of the approach used. There 
are some more or less rigorous theorems in mathematical analysis about dependence of ODE's 
solutions on parameter (see e.g. Korn & Korn (1968), sect. 10.2-7c, and Kamke (1959), 
sect. 2.5., SK78 and references therein) and this note is not the right place to discuss them, 
see e.g. SK78. In our case, we loosely reduce these theorems to assertion that, if for n = 
and n < 1 the solutions of ODE exist, then the formal expansion of ODE and its solution in 
the above mentioned way converge to the correct solution (actually, this is rewording of the 
Poincare-Lyapunov theorem, see SK78 and references therein). In our case of LEE we "do 
know" that desired solution exists (and we may find it e.g. by the numerical integration of 
ODE) therefore we may conclude that the procedure should give (in principle - for infinite 
number of terms in eq. (5)!) correct solution. Does the particular truncated approximation 
give correct (at least in the numerical sense) solution - this is not guaranteed and should be 
checked each time separately. 

As it happens in our case, one meets no trouble with the perturbation method in the first 
approximation, to first order of n (see however sect. 7), while already in the second approx- 
imation there are some difficulties. 



4. First approximation 



First approximation (zero'th approximation coincides with LEE of index n = and 
needs no additional solving) was solved in SK78 and we present here results with some 
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additional formulas needed for the second approximation. 

We have from eq. (6) (I remind that the shorthand s = \/6 is used): 

y (x) = 1 - - x 2 ; y (s) = 0; y' (s) = -- s; y' \s) = ~. 
Using eqs. (6, 12) we get from eq. (7): 

y[{x) = -\ [ t 2 \n(y (t))dt, 
x Jo 



or 



,..42 1 . . N 2s , 
VA X ) = - + -R x m (l/o) 5 -arctanh(- / 

x 9 x x 2 s 



We need also value of y[ at x = s: 



Vi(s) = -s(4-31n2). 



Using eq. (14) we get from eq. (8): 



yi{x) = / y'i(t))dt, 



or 



5 5 3/ 

yi(x) = — -4+(2 + yo) ln(y ) + 4 - arctanh (-), 

J. o 37 5 

j/i(s) = 41n2- \. 



(12) 

(13) 
(14) 

(15) 

(16) 

(17) 
(18) 



5. Second approximation 



We get from eq. (8): 



1 



y 2 (aO = — / * 

x Jo 



2/i (*) , 1 ,_2 



2/o(t) 2 



+ J in 2 (yo(0) 



(19) 



Using yo(^) from eq. (12) and yi(t) from eq. (17) we can integrate eq. (19) analytically 
however expression is cumbersome and we do not present it here. We only mention that 
y' 2 (x) has a log- singularity at x = s: 



Vmy 2 (x) = -\ /6 2(1 ^ /s) dx= S - yi (s) ln(l-^). 



(20) 
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Using solution y' 2 (x) we get further: 

V2{x)= I y 2 {t)dt. (21) 



Again analytical solution is very cumbersome and we present only the value of y2(x) at x = s: 

y 2 (s) = - (413 - 21 7T 2 - 402 ln(2) + 144 ln 2 (2)) . (22) 
9 



6. Series solution 

In SK77 the series solution was given for LEE of arbitrary n: 

oo 

y(x) = l + Y, a k x 2k , (23) 
k=i 

with recursion relation between coefficients a^: 

^ m 

a m +i = —, ' T\7r\ r^T y^(in + i - m)(m - i + 1) (3 + 2 (m - i)) aia m+1 _ h (24) 

i=l 

where m > 1, and a\ = —1/6. This recursion relation is valid for any n. The important 
problem of the convergence radius r of the series (23, 24) can not be solved in the general 
case. 

We consider here the interesting case of small values of n. In the linear approximation, we 
have for the first several coefficients: 

1 n 5n 70n 3150n . . 

ai = ~3!' a2= 5!' ° 3 = 3T7!' ° 4 = 9^9!' ° 5 = 45^1! " (25) 

We note that starting with a 2 all coefficients are positive. To derive a general formula for 
coefficients (in the same approximation, to the first order of n), we retain in the sum in r.s. 
of the recursion relation (24) only first (with i — 1) and last (with % = m) terms, and we get 
the following relation: 

(m-l)(2m+l) n 1 n 

a m+ i = — — — — —^a m , m>2, a x = --, a 2 = — 26 

6(m+l) (2m + 3) 6 5! 

From this we get for the convergence radius of the series solution of the LEE for quasi- 
incompressible planet: 

r= i im n^ =s . (27) 

m^oo y a m+ i 
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7. Density and pressure distribution 

Besides solving LEE (which is of its own interest and import) it's interesting to look 
for distribution of some physical parameters, namely density and pressure over the model's 
volume. We mention that under the gravitationally equilibrium condition, the pressure dis- 
tribution is more strictly defined: the pressure should be a continuous function of radius 
(while density can have discontinuity - as in the first order phase transition) and the pres- 
sure should be monotonically decreasing. 

Consider the spherical polytrope of index n, with central pressure P c , and central density p c . 
Transition from the "physical" current radius (distance from the center), r, to the "dimen- 
sional" radius x (which appears in all equations above) is as follows: 



a ■ x, a 



<" +1 > P <l 1/2 , (28) 



4%G p\ 



where G is the gravitational constant. 

Pressure P(r) and density p (r) at distance r are given by: 

P (r)/P c = y 1+n (x), p (r)/p c = y n (x), (29) 

where y (x) is the function which appears in all equations above. 

In the incompressible planet (polytropic index n — 0), the density p is constant over the 
volume, while the pressure decreases as: 

P(r) = P c -y (x) = P c ■ (1 - (x/s) 2 ). (30) 

In the quasi- incompressible planet, in the first approximation (linear in n) for density as a 
function of radius we have: 

P (r)/Pc = y n (x) = 1 + n In (y (x)) = 1 + n In (1 - (x/s) 2 ), (31) 

and at x = s (inside the planet close to its surface) there is a log-singularity of density. This 
is the reducible singularity, that is by integrating over volume we get correct values of radius, 
mass, and central-to-mean density ratio. Still this singularity in the density distribution is 
the unpleasant and non-physical characteristic of the perturbation method. Sure, in reality 
there is no density singularity inside the quasi-incompressible planet, and we may avoid it 
in our calculation keeping n in power of y(x) and rewrite eq. (31) as follows: 

P (r)/Pc = y n (x) = [y (x) +n yi (x)] n , (32) 

or even more accurately: 

p(r)/ Pc = y n (x) = [yo(x)+n yi (x)+n 2 y 2 (x)} n . (33) 
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As to the pressure distribution, there is no difficulty because in the linear approximation we 
have: 

P(r)/P c = y n+ \x) = y (x) + n [ yi (x) + y (x) In (y (x))} , (34) 

and there is no singularity, in the pressure distribution, inside the model (while there is a 
log- singularity in derivative dP/ dr which is proportional to the density - in accordance with 
eq. (31)). 

Again, as for the density distribution, we may write more accurately: 

P{r)/P c = [yo(x)+n yi (x)+n 2 y 2 (x)} n+1 , (35) 
avoiding the singularity in the derivative dP/dr. 

8. Global parameters 

Now we are ready to find some global parameters of the quasi-incompressible planet. 



8.1. Radius 



According to eq. (28), the total radius R of the polytrope is 

R = a ■ A, (36) 

where the first zero A is such that: 

y (X) = y (X) +n Vl (X) + n 2 y 2 (X) = 0, (37) 
X = s + nSi+n 2 S 2 , (38) 

with Si, S 2 still to be found. 

In the spirit of the perturbation method, we should keep, in eq. (37), the argument X in 
(A) to the second order of n, in yi(X), to the first order of n, and in y 2 (A), only A = s. 
Expanding all functions to the second order of n, combining terms with the same power of 
n, equalizing coefficients of n and n 2 to zero, we get system of two equations for Si and S 2 , 
and solving these equations we get finally: 
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§ 2 = — (1379 - 84 7T 2 - 888 In 2 + 144 In 2 2) . (41) 
From eqs. (38, 39, 41) we have numerically: 

X = 2.44948974278 + 0.537975784794 n + 0.12328309 n 2 . (42) 

Note that X (and hence radius R) is the only parameter of the quasi-incompressible planet 
which can be calculated to the second order of n. Due to the singularity of y' 2 (s), the second 
approximation is not possible for other global parameters of the quasi-incompressible planet. 



8.1.1. Comparison with the numerical values 

For n = 1/2, numerically X nurn = 2.752698, analytically, from eq. (42) X = 2.749298 
with relative difference of about -.12%. 

For n = 1/10, numerically X num = 2.504545, analytically X = 2.50452015 with relative 
difference of .01%. 

For n = 1/1000, numerically X num = 2.45004, analytically X = 2.45002978 with difference 
in the last digit of X num . 

We mention that the cases of very small values of n are more easier (and more accurately) 
described by analytical formulas than by numerical calculations. 

Besides, no numerical calculation can give analytical formula for X{n) even in the linear 
approximation in n. 



8.2. Mass 

The running mass inside the polytropic sphere of radius r is: 

m(r) = 4 7r p c a 3 [—x 2 y (x)] , 
and the total mass of the polytrope is: 

M = 4vr p c a 3 fi, // = -X 2 y'{X). 



(43) 



(44) 



In order to calculate /i(X) to the second order of n, we need y' 2 (s) which is infinite (see eq. 
(20)), so we restrict ourselves by the linear approximation, already considered in SK78. By 
the same way as in section (8.1), we get: 



/i = 2 s + 6 n 



S 1 -y' 1 ( s ) = 2s [1 - (37/6-8 In 2) n] = 2s(l - 0.62148922 n). (45) 
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8.3. Central-to-mean density ratio 

For the polytrope, from eqs. (36, 44) we have for the central-to-mean density ratio: 

Pm 3/i 

and though we calculated X in the second approximation, however the parameter /j, was 
calculated only in the first approximation, therefore we can calculate the central-to-mean 
density ratio for the quasi-incompressible planet only in the first approximation, already 
considered in SK78: 



— = 1 + - - 2 In 2 n = 1 + 1.280372 n. (47) 

Pm V 3 



8.4. Moment of inertia 

The central moment of inertia of the spherical polytrope is: 

nM i rX 

I=j r 2 dm = ki M R 2 , fc 7 = — - / y n (x)x A dx, (48) 
Jo pX Jo 

and to the first order of n we get: 

5 ~ 25' (49) 



8.5. Milne integral 



In the theory of polytropes, there is Milne's relation (Milne 1929) which also appeared 
in the "refined" theory of rotating polytropes ((Chandrasekhar & Lebovitz 1962), eq. (87)): 



MI 



= / 

Jo 



y 1+n {x) x 2 dx = 



n 



n 



(50) 



Using our formulas we find for both sides of the equation (50), for the quasi-incompressible 
planet (n « 1): 

MI = - s + — s (120 In 2 — 79) n. (51) 
5 25 

We can not calculate MI in the second approximation because this requires the y' 2 (s) (see 
r.s. of eq. (50)) which diverges (see eq. (20)). 

By the way, the case of n = 5 in eq. (50) is interesting as well. Assuming (5 — nO as a 



-10- 



small parameter, expanding both sides of eq. (50) in series we find the dependence of X(n), 
already found in SK78: 

X(n) = ^5— *— , 0<5-ra«l. (52) 
7r 5 — n 



8.6. One interesting relation 

In ( Seidov & Skvirsky 2000), the interesting relation is introduced between the grav- 
itational potential energy, W, the central potential, U c , and the mass of the celestial body. 
We confine ourselves here by the case of spherical polytropes. 

From the theory of polytropes, the central gravitational potential, U c , the central pressure 
P c , and the central density, p c , are (see (Chandrasekhar 1957), p. 100, eq. (85); p. 99, eq. 
(80,81); p. 78, eq. (99)): 

£/ c = (l + n)- + ^, (53) 
p — I ^1 (54) 

1 X X M 

Pc 3[-y'(X)] Pm [-y'(X)]4nRf 

Combining eqs. (53, 54, 55), we get for the central gravitational potential of the spherical 
polytrope: 

Potential energy of polytrope (we loosely take positive sign) is ((Chandrasekhar 1957), p. 
101, eq. (90)): 

3 G M 2 . . 

H '~— < 57 > 

Now we have, for polytropes, the VFf/M-ratio introduced in ( Seidov & Skvirsky 2000): 

W 3 1 

WUM = — — = — —— . (58) 

U C M 5 - n 1 + X/p K ' 

Note in brackets that WU ^M-ratio is remarkably constant (=2/5) for the rotating triaxial ho- 
mogeneous ellipsoids (see ( Seidov & Skvirsky 2000)). In our case of the quasi-incompressible 
planet, using eqs. (38, 39, 45, 58) we have to the first order of n: 



WUM = - 

5 



22 

1 -n{ 2 In 2) 



(59) 
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8.7. Parameter uj 



In (Christensen-Dalsgaard & Mullan 1994), the parameter uj for the polytrope was 
introduced (see their eq. (A14)): 



1 



UJ 



, i) 

4vr(l + n) \XJ 



l+l/n 



-i/'PO 



l-l/n 



(60) 



with a note that at n = this expression is undefined, and from numerical calculations, at 
n = 0, uj « 0.033175949. This is a good chance to demonstrate usefulness of the perturbation 
method by SK: using the formulas of the previous sections and omitting some algebra we 
write down the final expression for uj to the first order of n: 



UJ 



2tt 



exp(-8/3) < 1 -n 



1 + 



1 

2 \ 3 



In 6 



(61) 



or numerically 



uj = 0.033175904175 (1. - 1.382731302 n), 



and we note again that the analytical expression helps to check the numerical calculations 
(and the last two digits of the "numerical" value of uj at n = differ from the "analytical" 
value) . 



9. "Numerical" perturbation method 

Though being slightly out-of-topic, we mention that the perturbation method can be 
used in its numerical modification as well. 

In the case of polytropes, it is particularly interesting to do this for the case n = 3. We 
briefly describe the method and give the short results in the first approximation. We take 
in the r.s. of eq. (1), n — 3 — 5 with < 5 <C 3, then we have 

(yo - 5 yif~ 5 = yl~yl (3 y x + y In y ) 5. (62) 

Here yo is solution of LEE of index n = 3 and y\ is the "perturbation" function. Also, 
because S > 0, zero X of function yo — 5 y± is less than Xo, zero of function yo. Due to an 
inequality X < X there is no problems with behavior of functions near boundary as in the 
quasi-incompressible case. As in the text above, we have system of two ODE's (each of them 
of the second order) for functions y and y±. We can solve this system numerically and find 
all relevant functions yo, yi, y' , and y[ and particularly their values at x = X . Then we can 
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get the formulas similar to ones in section (8). 

We present the results of the numerical calculations: 



X = 6.8968498, y' {X ) = -0.0424297317, 



(63) 



yi(X ) = 0.16547670449, yi (X ) = -0.00616735. 



(64) 



From this we have for the dimensionless radius of polytrope with index n close to 3: 



Numerically, atn = 3 — 0.1, X — 6.526374, while from the linear perturbation method (65): 
X = 6.5068. 

Also, numerically, ai n — 3 + 0.1, X — 7.308484, while from the linear perturbation method 
(65): X = 7.2868. 

Deviations of numerical values in two cases correspond (correctly) to the positive second 
derivative of function X{n). 



In this note we present the exact analytical solutions for the internal structure and 
global parameters of the quasi-incompressible planet modeled as the polytrope of small 
index < n C 1. The perturbation method used here is not rigorously justified by means 
of the theory of differential equations and there are some problems about application of the 
method in the interval of argument where the perturbation function is of the same order 
or even larger than the initial non-perturbed function. The problem of justification is here 
similar to the problem of rotationally distorted polytropes which was already discussed some 
in the astrophysical literature (see e.g. Chandrasekhar & Lebovitz (1962)). Still validity of 
any method in applicational sciences (as astrophysics is such relative to mathematics) may 
be compared by numerical calculation and (astro)-physical "common sense" and I only hope 
that this note may trigger some discussion as well. 



This paper is to be devoted to the memory of my dear friend and PhD student late Rafael 
Khejrullaevich Kuzakhmedov (=Rafik) who having been a very strong theor. physicist flatly 
rejected to admit any significance in our two minor notes (SK77, SK78) and had never got 
even PhD degree. Only my vision of astrophysics and Rafik's incredible ability of what is 



X = 6.896848 + (3 - n) y 1 (X )/y' (X ) = 6.896848 - 3.90002 (3 - n). 



(65) 



10. Summary 
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